"""
# LidDrivenCavity.jl

A Julia package for solving lid-driven cavity flow problems using the stream function-vorticity formulation
on L-shaped domains. Features integration with the GCR.jl framework for high-performance linear algebra
and optional PETSc backend support for parallel computing.

## Key Features

- Stream function-vorticity formulation for 2D and 3D incompressible Navier-Stokes equations
- L-shaped domain handling with complex boundary conditions
- Integration with GCR.jl for advanced iterative linear solvers
- Optional PETSc backend for parallel computing and scalability
- Automatic relaxation parameter adjustment for numerical stability
- Comprehensive convergence analysis and visualization tools

## Solver Types

- **2D Solver**: Classic stream function ψ and scalar vorticity ω formulation
- **3D Solver**: Vector stream function (ψₓ, ψᵧ, ψᵤ) and vector vorticity (ωₓ, ωᵧ, ωᵤ) formulation

## Linear Solver Backends

- **Pure Julia**: Native GCR implementation for moderate-sized problems
- **PETSc**: Distributed parallel solvers for large-scale HPC applications
- **GPU**: CUDA/ROCm acceleration for supported hardware

"""
module LidDrivenCavity

# Core Julia dependencies
using LinearAlgebra
using SparseArrays
using Printf
using IterativeSolvers
using Statistics
using Random

# Enhanced GCR.jl availability detection with multiple fallback paths
const HAS_GCR, GCR_MODULE = begin
    local has_gcr = false
    local gcr_mod = nothing
    
    try
        # Try multiple paths for GCR.jl
        gcr_paths = [
            "/home/linden/code/work/Helmholtz/gcr-nccl/GCR-Julia/src",
            "/home/linden/code/work/Helmholtz/gcr-nccl/GCR-Julia",
            joinpath(dirname(@__DIR__), "..", "..", "work", "Helmholtz", "gcr-nccl", "GCR-Julia", "src"),
            joinpath(dirname(@__DIR__), "..", "..", "work", "Helmholtz", "gcr-nccl", "GCR-Julia")
        ]
        
        # Add paths that exist to LOAD_PATH
        for path in gcr_paths
            if isdir(path) && !(path in LOAD_PATH)
                push!(LOAD_PATH, path)
            end
        end
        
        # Try to import GCR module
        gcr_mod = Base.require(Main, :GCR)
        has_gcr = true
        @info "Successfully loaded GCR.jl framework from dynamic path detection"
        
    catch e1
        # Fallback: try direct using statement
        try
            using GCR
            gcr_mod = GCR
            has_gcr = true
            @info "Successfully loaded GCR.jl framework via direct import"
        catch e2
            @warn "Could not load GCR.jl via any method: Primary: $e1, Fallback: $e2. Using fallback solvers."
            has_gcr = false
            gcr_mod = nothing
        end
    end
    
    (has_gcr, gcr_mod)
end

# Export main solver interfaces
export solve_lid_driven_2d, solve_lid_driven_3d
export LidDrivenOptions, LidDrivenResult
export LShapeDomain, create_lshape_domain

# Export geometry and grid utilities
export build_grid, fluid_mask, boundary_mask, interior_mask
export apply_boundary_conditions!

# Export numerical method components
export assemble_poisson_matrix, finite_difference_operators
export vorticity_transport!, boundary_vorticity_update!
export under_relax!, auto_adjust_relaxation

# Export analysis and diagnostics
export convergence_analysis, plot_solution, validate_solution
export performance_metrics, compare_solvers

# Export linear algebra enhancements
export GCRPoissonSolver, FallbackPoissonSolver, AdaptivePoissonSolver
export PETScPoissonSolver, create_petsc_poisson_solver
export assemble_poisson_matrix_efficient
export comprehensive_solver_benchmark, analyze_benchmark_results
export MemoryTracker, get_memory_usage, update_memory_tracker!

# Define package version and constants
const VERSION = v"0.1.0"
const REYNOLDS_DEFAULT = 100.0
const TOLERANCE_DEFAULT = 1e-5

"""
    LidDrivenOptions

Configuration options for lid-driven cavity solvers.

# Fields
- `n::Int`: Number of interior grid intervals per axis (default: 32)
- `Re::Float64`: Reynolds number (default: 100.0)  
- `dt::Union{Float64,Nothing}`: Time step (auto-computed if nothing)
- `max_steps::Int`: Maximum iterations (default: 10000)
- `tol::Float64`: Convergence tolerance (default: 1e-5)
- `α_ψ::Float64`: Stream function relaxation factor (default: auto)
- `α_ω::Float64`: Vorticity relaxation factor (default: auto)
- `solver::Symbol`: Linear solver backend (`:julia`, `:petsc`, `:gpu`)
- `preconditioner::Symbol`: Preconditioner type (`:none`, `:diagonal`, `:ilu`)
- `verbose::Bool`: Enable detailed output (default: false)
- `save_history::Bool`: Save convergence history (default: true)
- `U_lid::Float64`: Lid velocity magnitude (default: 1.0)
"""
Base.@kwdef struct LidDrivenOptions
    n::Int = 32
    Re::Float64 = REYNOLDS_DEFAULT  
    dt::Union{Float64,Nothing} = nothing
    max_steps::Int = 10000
    tol::Float64 = TOLERANCE_DEFAULT
    α_ψ::Union{Float64,Nothing} = nothing  # Auto-adjust if nothing
    α_ω::Union{Float64,Nothing} = nothing  # Auto-adjust if nothing  
    solver::Symbol = :julia
    preconditioner::Symbol = :diagonal
    verbose::Bool = false
    save_history::Bool = true
    U_lid::Float64 = 1.0
end

"""
    LidDrivenResult

Results from lid-driven cavity solver.

# Fields
- `ψ::Array`: Stream function field(s)
- `ω::Array`: Vorticity field(s) 
- `u::Array`: x-velocity field
- `v::Array`: y-velocity field  
- `w::Union{Array,Nothing}`: z-velocity field (3D only)
- `x::Vector{Float64}`: x-coordinates
- `y::Vector{Float64}`: y-coordinates
- `z::Union{Vector{Float64},Nothing}`: z-coordinates (3D only)
- `converged::Bool`: Whether solver converged
- `iterations::Int`: Number of iterations performed
- `residual_norm::Float64`: Final residual norm
- `solve_time::Float64`: Total solve time (seconds)
- `convergence_history::Vector{Float64}`: Residual history per iteration
- `options::LidDrivenOptions`: Solver options used
"""
struct LidDrivenResult
    ψ::Array  
    ω::Array
    u::Array
    v::Array
    w::Union{Array,Nothing}
    x::Vector{Float64}
    y::Vector{Float64}
    z::Union{Vector{Float64},Nothing}
    converged::Bool
    iterations::Int
    residual_norm::Float64
    solve_time::Float64
    convergence_history::Vector{Float64}
    options::LidDrivenOptions
end

"""
    LShapeDomain{D}

Represents an L-shaped computational domain in D dimensions.

The L-shape is formed by removing the upper-right quadrant from a rectangular domain,
i.e., Ω = [-1,1]^D \\\\ [0,1]^D.

# Fields
- `bounds::NTuple{D,Tuple{Float64,Float64}}`: Domain bounds for each dimension  
- `removed_bounds::NTuple{D,Tuple{Float64,Float64}}`: Bounds of removed region
- `fluid_mask::Array{Bool,D}`: Mask indicating fluid nodes
- `boundary_mask::Array{Bool,D}`: Mask indicating boundary nodes
- `interior_mask::Array{Bool,D}`: Mask indicating interior nodes
"""
struct LShapeDomain{D}
    bounds::NTuple{D,Tuple{Float64,Float64}}
    removed_bounds::NTuple{D,Tuple{Float64,Float64}}  
    fluid_mask::Array{Bool,D}
    boundary_mask::Array{Bool,D}
    interior_mask::Array{Bool,D}
end

# Include core implementations
include("geometry/lshape_domain.jl")       # L-shape domain geometry
include("geometry/grid_generation.jl")     # Grid generation utilities  
include("geometry/boundary_conditions.jl") # Boundary condition handling

include("numerics/finite_differences.jl")  # Discrete operators
include("numerics/time_integration.jl")    # Time stepping schemes
include("numerics/relaxation.jl")          # Under-relaxation methods

include("linear_algebra/poisson_assembly.jl")  # Poisson matrix assembly
include("linear_algebra/gcr_integration.jl")   # GCR.jl integration layer
include("linear_algebra/petsc_backend.jl")     # PETSc backend integration
include("linear_algebra/performance_monitoring.jl")  # Performance and memory monitoring

include("solvers/navier_stokes_2d.jl")     # 2D stream function-vorticity solver
include("solvers/navier_stokes_3d.jl")     # 3D vector stream function solver

include("analysis/convergence_analysis.jl") # Convergence diagnostics
include("analysis/validation.jl")          # Solution validation
include("analysis/visualization.jl")       # Plotting and visualization

"""
    solve_lid_driven_2d(options::LidDrivenOptions=LidDrivenOptions()) -> LidDrivenResult

Solve the 2D lid-driven cavity problem on an L-shaped domain using stream function-vorticity formulation.

The problem solves the coupled system:
- Poisson equation: ∇²ψ = -ω (stream function)  
- Vorticity transport: ∂ω/∂t + u·∇ω = (1/Re)∇²ω (pseudo-time stepping to steady state)

# Arguments
- `options::LidDrivenOptions`: Solver configuration options

# Returns
- `LidDrivenResult`: Complete solution with fields, diagnostics, and metadata

# Example
```julia
# Basic usage with default parameters
result = solve_lid_driven_2d()

# Custom configuration  
opts = LidDrivenOptions(n=64, Re=1000.0, solver=:petsc, verbose=true)
result = solve_lid_driven_2d(opts)

# Access solution fields
ψ, ω = result.ψ, result.ω
u, v = result.u, result.v
println("Converged: \$(result.converged) in \$(result.iterations) iterations")
```
"""
function solve_lid_driven_2d(options::LidDrivenOptions=LidDrivenOptions())
    return solve_lid_driven_2d_impl(options)
end

"""
    solve_lid_driven_3d(options::LidDrivenOptions=LidDrivenOptions()) -> LidDrivenResult

Solve the 3D lid-driven cavity problem on an L-shaped domain using vector stream function-vorticity formulation.

The problem solves the coupled system:
- Poisson equations: ∇²ψₖ = -ωₖ for k ∈ {x,y,z} (vector stream function)
- Vorticity transport: ∂ω/∂t + u·∇ω - ω·∇u = (1/Re)∇²ω (includes stretching term)

# Arguments  
- `options::LidDrivenOptions`: Solver configuration options

# Returns
- `LidDrivenResult`: Complete solution with vector fields and 3D diagnostics

# Example
```julia
# 3D simulation (computationally intensive)
opts = LidDrivenOptions(n=16, Re=100.0, max_steps=500)
result = solve_lid_driven_3d(opts)

# Access 3D vector fields
ψx, ψy, ψz = result.ψ[:,:,:,1], result.ψ[:,:,:,2], result.ψ[:,:,:,3]  
ωx, ωy, ωz = result.ω[:,:,:,1], result.ω[:,:,:,2], result.ω[:,:,:,3]
u, v, w = result.u, result.v, result.w
```
"""  
function solve_lid_driven_3d(options::LidDrivenOptions=LidDrivenOptions())
    return solve_lid_driven_3d_impl(options)
end

"""
    create_lshape_domain(n::Int, ::Val{D}=Val(2)) -> LShapeDomain{D}

Create an L-shaped domain with n grid intervals per axis.

# Arguments
- `n::Int`: Number of interior grid intervals per axis
- `Val{D}`: Dimension (2 or 3)

# Returns  
- `LShapeDomain{D}`: L-shaped domain with precomputed masks
"""
function create_lshape_domain(n::Int, ::Val{D}=Val(2)) where D
    return create_lshape_domain_impl(n, Val(D))
end

# Package information functions
"""
    lid_driven_info()

Display information about the LidDrivenCavity.jl package installation and capabilities.
"""
function lid_driven_info()
    println("LidDrivenCavity.jl v$(VERSION)")
    println("High-Performance Lid-Driven Cavity Flow Solver")  
    println()
    
    println("Supported Features:")
    println("  ✓ 2D and 3D L-shaped domains")
    println("  ✓ Stream function-vorticity formulation")
    println("  ✓ Automatic relaxation parameter tuning")
    println("  ✓ Multiple linear solver backends")
    
    if HAS_GCR
        println("  ✓ GCR.jl framework integration")
        println("  ✓ Advanced iterative solvers (GCR, CA-GCR, BCA-GMRES)")
        try
            # Check for PETSc extension
            println("  ✓ PETSc parallel computing support")  
        catch
            println("  ○ PETSc support (requires PETScExt.jl)")
        end
    else
        println("  ○ GCR.jl integration (not available)")
        println("  ○ Advanced solver algorithms (requires GCR.jl)")
    end
    
    println()
    println("Example usage:")
    println("  julia> using LidDrivenCavity")
    println("  julia> result = solve_lid_driven_2d()")
    println("  julia> plot_solution(result)")
end

# Module initialization  
function __init__()
    if HAS_GCR
        @info "LidDrivenCavity.jl initialized with GCR.jl framework"
    else
        @warn "LidDrivenCavity.jl initialized without GCR.jl (limited functionality)"  
    end
end

end # module LidDrivenCavity